New results for electromagnetic quasinormal modes of black holes 
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The differential equations governing the late-time ring-down of the perturbations of the Kerr 
metric, the Teukolsky Angular Equation and the Teukolsky Radial Equation, can be solved analyt- 
ically in terms of confluent Heun functions. In this article, for the first time, we use those exact 
solutions to obtain the electromagnetic (EM) quasinormal spectra of the Kerr black hole . This is 
done by imposing the appropriate boundary conditions on the solutions and solving numerically the 
so-obtained two-dimensional transcendental system. 

The EM QNM spectra are compared with already published results, evaluated through the con- 
tinued fraction method. The comparison shows that the modes with lower n coincide for both 
methods, while those with higher n may demonstrate significant deviations. To study those devi- 
ations, we employ the e-method which introduces small variations in the argument of the complex 
radial variable. Using the e-method one can move in the complex r— plane the branch cuts of the 
solutions of the radial equation and examine the dependence of the spectrum on their position. 

For different values of e, one can obtain both the frequencies evaluated through the well-established 
continued fraction method or somewhat different spectra calculated here for the first time. This 
result raises the question what spectrum should be compared with the observational data and 
why. This choice should come from better understanding of the physics of the problem and it may 
become particularly important considering the recent interest in the spectra of the electromagnetic 
counterparts of events producing gravitational waves. 



QUASI-NORMAL MODES OF BLACK HOLES 

During the long history of the study of the quasi- 
normal modes (QNMs) of a black hole (BH) ([l-oO]), 
the case of electromagnetic (EM) perturbations has been 
often ignored because of their perceived irrelevance to 
the problem of finding gravitational waves. The reasons 
behind this are that: 

1. The expected luminosity of the gravitational output 
in the most often studied process of BH binary merger is 
much bigger than the electromagnetic one ([■Vi]). 

2. The EM waves strongly interact with the surround- 
ing medium. This may lead to essential deviation of the 
observed spectrum from the one expected from the no- 
hair theorem and makes the object's fingerprint harder 
to detect. 

3. Most importantly, because of this interaction, the 
electromagnetic perturbations are strongly absorbed by 
the interstellar medium, thus making the detection of the 
signal almost impossible at the predicted low frequencies 
for the electromagnetic QNMs. 

On the other hand, the gravitational waves (GW) in- 
teract very weakly with matter and thus they can be 
detected at big distances, without getting absorbed or 
scattered, i.e. without obscuring the signature of the 
body that emitted them. It is, therefore, reasonable to 
expect that the GW should be much better suited for 
studying the central engine of astrophysical events, such 
as gamma-ray bursts (GRBs), while the EM waves should 
be seriously influenced by its eirvironment. 

For now, however, there are no gravitational waves de- 
tected. Although both LIGO and VIRGO detectors al- 



ready work at design sensitivity, both detectors still fail 
to "see" gravitational waves ([32-38]). Particularly puz- 
zling is the lack of GW detection from short GRBs ( 
[39, 40]) whose progenitors are expected to emit GWs in 
the range of sensitivity of the detectors. 

The simplest explanation of those negative results may 
be a new mechanism of generation of short GRBs which 
in good approximation preserves the spherical symme- 
try of the process in the central engine and admits only 
significant dipole radiation. As a result, no significant 
gravitational waves will be generated during the short 
GRBs, since the gravitational waves have a quadrupole 
character. A similar situation is observed in the long 
GRBs. Such a new hypothesis for short GRBs is sup- 
ported by the strong observational indications that both 
types of GRBs may have the same nature and differ only 
in their time scales ([41, 42]). If so, we may expect that 
most of the energy release from GRBs is in the form of 
electromagnetic radiation. 

A more traditional point of view is a physical process 
which yields both electromagnetic and gravitational ra- 
diation from GRBs. Actually, the ratio of the energy 
release in the form of electromagnetic waves and in the 
form of gravitational waves is still an open problem. Its 
solution strongly depends on the details of the hypothet- 
ical mechanism of GRBs which is still far from being well 
established. For example, the expected (but not yet ob- 
served) energy output in the case of GWs from a BH 
merger is ^ lO^^erg, which coincides (up to a factor due 
to coUimation) with the observed electromagnetic energy 
output of GRBs. 

While hopes are laid on the Advanced LIGO and Ad- 
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vanced VIRGO which should start operating in the next 
years, this situation offers a good motivation for opti- 
mizing the GW search strategy and understanding bet- 
ter the physics of the GW sources. Particularly, this 
points to the advantages of studying the EM counter- 
part of the GW emission, which can help the localiza- 
tion of the source (improving on the big error box of the 
GW detector) and also it may give additional clues to 
the physics of the event ([31, 43]). Numerical simula- 
tions already explore the detectability of the EM coun- 
terpart in different cases (for the case of supermassive 
black holes mergers, for example, see [45, 46], for neu- 
tron stars mergers ]43, 44]). The first results of LIGO 
and VIRGO searches using such a multimessenger ap- 
proach were also published [47]. The idea behind those 
searches (for details, see [47, 4.S]) is to use the GRBs 
as the EM counterparts of the GW, since the suspected 
GRB progenitors (collapsars for long GRBs and binary 
system mergers of neutron stars and/or black holes for 
the short GRBs) should emit GWs as well as EM and 
there is already a well working mechanism for observing 
the extremely EM luminous GRBs. Although the theo- 
retical results from the multimessenger approach are still 
being analyzed, the intensive activity in this field shows 
that the EM counterpart of the GW emission can both 
facilitate and improve the information obtained from the 
GW observations. 

The discrete spectrum of complex frequencies called 
QNMs describes only the linearized perturbations of the 
metric. Hence, they cannot describe completely the dy- 
namics of the process during the early, highly intensive 
period of those events when the linearized theory is not 
applicable. On the other hand, it is known from full nu- 
merical simulations that it is the QNMs which dominate 
the late-time evolution of the object response to pertur- 
bations ]19]. Thus, from an observational point of view, 
the QNM are important, since we may observe only the 
tails from the corresponding events, being far from them. 
This conclusion is supported by the recent numerical ob- 
servation of two lowest gravitational modes of QNMs in 
the spectrum of the signal obtained from the full 3d gen- 
eral relativistic head-on collision of non-spinning BH (see 
]49], and for further information [^>()]). This result is not 
isolated - there are number of works in which the QNMs 
approximate well the signal of full 3d general relativistic 
simulations of mergers (for example [''i l-^ri] and also the 
pioneer works discussed in ]-'■">]). This clearly implies that 
studying QNMs can bring new insights to the physics in 
the processes which include strong-field regime. 

Those numerical results also point to another pos- 
sible use of the QNMs in astrophysical observations. 
The QNMs correspond to particular boundary conditions 
characteristic for the object in question, and since in the 
case of BH the no-hair theorem states that they should 
depend only on the parameters of the metric (the mass 
AI and the rotation a for the case of Kerr BH), mea- 



suring those frequencies observationally can be used to 
test the nature of the object - a black hole or other com- 
pact massive objects like super-spinars (naked singulari- 
ties), neutron stars, black hole mimickers etc. [r)(i"()l]. It 
also can constrain additionally the no-hair theorem which 
was recently put into question in the case of black holes 
formed as a result from the collapse of rotating neutron 
stars ]()2]. An interesting possibility is to find a way to 
use the damping times of the EM quasi-normal modes for 
comparison with observations. While the frequencies are 
subject to interaction with the surrounding matter which 
can significantly change the spectra, their damping times 
should be much less prone to deviation. A suggestion for 
such use can be found in simulations of jet propagation, 
which imply that the short-scale variability of the light- 
curve should be due to the central engine and not to the 
interaction of the jet with the surrounding medium (see 
[63] and reference therein). 

One more important application of the QNM spectrum 
can be found in the study of the central engines of the 
GRBs, whose extreme luminosity(~ 10^^ — lO^^erg/ s) 
and peculiar time-variability cannot yet be fully ex- 
plained in the frames of current models. Even though 
numerical simulations proved to be capable of describ- 
ing some of the features of the GRBs light curves (for a 
recent review on GRBs, see [()4]), the biggest stumbling 
stone seems to be the lack of proper understanding of the 
central engine of the GRBs. 

Common ingredients of the existing GRB models in- 
clude a compact massive object (black hole or a milli- 
second magnetar) and extreme magnetic field lO^^G 
) which accelerate and coUimate the matter via differ- 
ent processes. Although those processes are still an open 
question for both theory and numerical simulations, the 
very central engine can be studied approximately by the 
linearized EM (and also GW if data is available) pertur- 
bations of the Kerr metric. When finding the electromag- 
netic QNMs, one does not care for the origin of the per- 
turbation, but only for its spin and the parameters of the 
compact massive object. In the idealized EM case, the 
perturbation is described by free EM waves in vacuum. 
While the astrophysical black holes are thought to be not 
charged, they are immersed into EM waves with differ- 
ent energy and origin. The black hole response to such 
EM perturbations in linear approximation will be then 
the QNM spectrum defined by the appropriate boundary 
conditions ^. 

Studying the so obtained electromagnetic spectrum 
can give important insights into the key parameters of 
the physics occurring during high-energy events as GRBs. 
In particular, the electromagnetic QNMs are subject to 



^ Other conditions more suitable for describing a primary jet were 
studied in [61]. 
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resonant amplification (the idea of tlie black hole bomb, 
[8, 28, 30]) and additionally, it is known from previous 
evaluations of the spectrum, that they exhibit very low 
damping in the limit a — > M. For the moment, there 
are no observations of the rotations of the GRB pro- 
genitors, but the theoretical expectations are that they 
should be highly rotating in order to produce jets with 
such luminosity and collimation. Available observation- 
ally measured rotations of astrophysical compact massive 
objects show that there are many cases of near extremal 
values thus studying the extremal limit could be rele- 
vant to such objects. For example, recent evaluations 
of the spin parameter of astrophysical black holes give 
for the spin parameter a = 0.63,0.90 and a = 0.89,0.99 
for AI = l,O.lA/0 for Sw J1644+57 and Sw J2058+05 
respectively (most probable values, see [65]), and also 
a > 0.98 for GRS 1915+105 (](iG]) and a = 0.989 for 
MCG-6-30-15 (]()?]). 

Moreover, because of the relatively good coverage of 
the GRBs observations, there is a great amount of data, 
in a wide energy range (from optical to GeV energies) 
and from different epochs of the bursts which can be used 
to test the eventual applicability of the QNM spectrum 
in the late-time epoch of the burst. It may be hard to 
extract EM QNM spectra from the existing crude GRB 
spectra since the basic EM QNM frequencies are very 
low (from a small part of Hz — for supermassive BH 
to several kHz - for stellar mass BH) and the intensity 
of the higher EM QNM may be very low. To the best 
of our knowledge such attempts haven't been made. A 
new space mission, which will additionally help the EM 
observations in the radio range - RadioAstron - will of- 
fer unprecedented resolution (up to 1/iarcsec) in a wide 
range of high frequencies (from 0.3GHz to 18-25GHz) ac- 
companied by continuum, polarized and spectral imaging 
(for details see RadioAstron website ^). One may hope 
to use this new mission for a more detailed study of the 
spectra of EM radiation from the compact objects but 
its sensitivity is also far from the area of the basic EM 
QNM. 

Theoretical calculations of the QNMs, however, are 
not simple. The linear perturbations of the rotating 
BHs are described by two second-order linear differen- 
tial equations: the Teukolsky radial equation (TRE) and 
the Teukolsky angular equation (TAE) on which specific 
boundary conditions should be imposed ([10, 15]). Until 
recently, solving those equations analytically was consid- 
ered impossible in terms of known functions, so approx- 
imations with simpler wave functions were used instead. 
The resulting system of spectral equations - a connected 
problem with two complex spectral parameters: the fre- 
quency Lu and the separation constant E - has been solved 



^ http:/ /www. asc.rssi.ru/radioastron 



using different methods ([20, 23, 25, 29]) with notably the 
most often used being the method of the continued frac- 
tion. This method was adapted by Leaver from the prob- 
lem of the hydrogen molecule ion in quantum mechanics 
]1(), 17]. While being successful in obtaining the QNMs 
spectra. Leaver's method has the disadvantage of not be- 
ing directly connected with the physics of the problem, 
thus making it harder to further explore the spectra ~ 
for example studying its dependence on the choice of the 
branch cuts of the exact solutions of the radial equation. 
In addition, one has some specific numerical problems in 
calculation of particular modes, for example, in calcula- 
tion of the 9*'' one in the gravitational case [18, 20, 25]. 

The analytical solutions of the TRE and the TAE can 
be written in terms of the confluent Heun function (for 
a 7^ M) as done for the first time in ]22, 26, 27, 68]. 
Those functions are the unique local Frobenius solu- 
tions of the second-order linear ordinary differential 
equation of the Fuchsian type [69-72] with 2 regular 
singularities {z = 0, 1) and one irregular (z = oo) 
(for details see [27]) and in the maple notation, they 
are denoted as: HeunC(Q!, /3, 7, (5, 77, 2) (normalized to 
HeunC(a, /3, 7, J, ry, 0) = 1). While the theory of the 
Heun functions is still far from being complete, they are 
implemented in the software package maple and despite 
the problems in that numerical realization (see the dis- 
cussion in [7.'^)]), the confluent Heun function was used 
successfully in our previous works [22, 61, 73, 74]. The 
advantage of using the analytical solutions is that one 
can impose the boundary conditions on them directly (see 
]22, 61]), and thus be able to control all the details of the 
physics of the problem. 

In a series of articles, we developed a method for solv- 
ing numerically two-dimensional systems featuring the 
Heun functions (the two-dimensional generalization of 
the Miiller method described in ]75]) and we used it suc- 
cessfully in the case of gravitational perturbation s = — 2 
of the Schwarzschild metric ]73]. The so obtained fre- 
quencies repeat with high precision the results already 
published by other authors. Additionally, we used the 
epsilon-method (see below) to study the branch cuts of 
the solutions, which are particularly important in the 
case of the 9*^* mode for s — —2. The latter, because 
of its very small real part, was often wrongly considered 
to represent the purely imaginary algebraically special 
mode. While the analysis of the potentials of the Regge- 
Wheeler equation (RWE) and the Zerilli equation (ZE) 
showed that there is a branch cut on the imaginary axis 
for this mode [7(), 77] which leads to its interesting prop- 
erties, this result is directly obtainable from the actual 
solutions of the RWE and ZE in terms of the confluent 
Heun functions. Furthermore, the stability of the solu- 
tions with respect to movement of the position of those 
branch cuts was studied in the whole interval of appli- 
cability of the method. Such a study is impossible with 
the continued fraction method, where the radial variable 
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does not explicitly enter the equations being solved and 
which cannot be used for purely imaginary frequencies 
([17], p. 8). If one looks at the equations used by this 
method in detail, it turns out that the angular equation 
[16] in the continued fraction method coincides with the 
the three-term recurrence defining the confluent Heun 
function, solution of the TAE, in the neighborhood of the 
two regular singular points, m = — 1, 1, where u = cos(0) 
([71] Eq. (1.9-1.10). The radial equation in the contin- 
ued fraction method, however, differs from the solution 
of the TRE in terms of the confluent Heun functions. 
This is because in Leaver's paper, the series from which 
the continued fraction are obtained are developed for the 
powers of (due to switching the places of the singu- 
lar points, see [17], p. 7), while the asymptotic three-term 
recurrence of the confluent Heun function at infinity is 
developed for ^J"^ . Note that in MAPLE, for r > r+, 
the evaluation of the confluent Heun functions at infin- 
ity is obtained by numerical integration from the second 
singularity r = r^. 

In this article, we continue the exploration of the ap- 
plication of the confluent Heun functions by studying the 
QNMs of the Kerr BH. Our results show that using the 
confluent Heun function, one can obtain the QNMs for 
a wide range of modes and rotational parameters, and 
that there is very good agreement between our results 
and those obtained within other methods. Using the e- 
method made it possible, for the first time to study the 
dependence of the so obtained frequencies on small de- 
viations in the phase condition and it is shown how this 
nontrivial dependence evolves with n and a. In this case, 
some of the modes are independent of e which should be 
expected since the frequencies should not depend on the 
radial variable. Other modes, however, depend critically 
on the value of e and they can differ seriously from the 
already published results. Additionally, details how the 
modes change in the interval of validity of the steepest 
descent method are presented. 



THE TEUKOLSKY ANGULAR EQUATION 

In Chandrasekhar's notation, the Teukolsky Master 
Equation ([>^]), for \s\ = 1 is separable under the sub- 



stitution ^ 



^S{e)R{r), where m = 0,±1,±2 
for integer spins and uj is the complex frequency. Be- 
cause of the choice of this form of XP, the sign of lu differs 
from the one Teukolsky used, and the stability condition, 
guaranteeing that the perturbations will damp with time, 
reads 3(w) > 0. 

The TAE for EM perturbations (s = —1) has 16 classes 
of exact solutions S{6) in terms of the confluent Heun 
functions (for full details see [27]). To fix the spec- 
trum approximately, one requires an additional regular- 
ity condition for the angular part of the perturbation, 
which means that if we choose one solution, 5*1 (9) regular 



around the one pole of the sphere {9 = 0) and another, 
82(9), which is regular around the other pole {9 = tt), 
then in order to ensure a simultaneous regularity, the 
Wronskian of the two solutions should become equal to 
zero, W[Si{9), S2{9)] = 0. This gives us one of the equa- 
tions for the two-dimensional system that needs to be 
solved to obtain the QNMs of the Kerr BH. 

In [27], there are four pairs of Wronskians, each pair 
being valid in a sector of the plane {s, m}. Ideally, using 
any of them should lead to the same spectrum. Numeri- 
cally, the results obtained with the different Wronskians 
coincide within 11-13 digits of precision. The Wronskians 
used to obtain the spectrum are: 



HeunC' (ai , /3i , 71 , (5i , r/i , (cos (7r/6) )^) 



HeunC(ai,/3i,7i,(5i,77i, (cos(7r/6)) ) 

HeunC'(a2,/32,72,<52,'f/2, (sin(7r/6))^) 

2 — \- P = ^ 

HeunC(a2,;52,72,<52,??2, (sin(7r/6)) ) 



(1) 



where the derivatives are with respect to z and the values 
of the parameters for the two confluent Heun functions 
for each m are as follows: 

For m = 0: ai ~ 4aa;,/3i = 1,71 = — l,(5i = 
4 aw, 771 = 1/2 — E — 2auj — a}uj^ and 

a.2 = -4aa;,/?2 = 1,72 = 1,(52 = -4 aw, 772 =,1/2- 
E + 2auj - d^u?, p = . . , 

' ' (sin(ir/6)) 

For m = 1: ai = — 4aw,/3i = 2,71 = 0,(5i = 
4 aio, i]i = 1 — E — 2 auj — (P'uP' and 

0L2 — —4 aw, /32 = 0, 72 = 2, ^2 = -4 aw, ,r]2 = I - E + 
2 aw — a^w^ and p — —4 aw 

For m ~ 2: a\ ^ — 4aw,/3i = 3, 71 = — l,(5i = 
4 aw, 771 = 5/2 — i? — 2 aw — a^w^ and 

a2 = -4 aw, /32 = 1, 72 = -3, ^2 = -4 aw, 772 = 5/2 - 
2 aw — c?ijj^ and p = 8 — 4aw. 
where we use 9 = ir/i (the QNMs should be independent 
of the choice of 9 in the spectral conditions) . 

These Wronskians differ from those in ]27], most no- 
tably by the presence of the term p. The reason for this is 
that they were constructed using different two solutions 
[5*1 (6*), 5*2(6')] of the TAE (note that the sign convention 
in this paper differs from the one in ]27]), each of which 
still being regular on one of the poles. That was done 
to improve the numerical convergence of the root-finding 
algorithm and to avoid maple's problems with the eval- 
uation of the confiuent Heun function and its derivative 
for certain values of the parameters. 



THE TEUKOLSKY RADIAL EQUATION 

The TRE differential equation is of the confluent Heun 
type, with r = r± regular singular points and r = 00 ~ 
irregular one. As it was noted in ]()1], the point r = 0, = 
7r/2 is not a singularity for this equation and, therefore. 
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it need not be considered when imposing the boundary 
conditions. The solutions of the TRE for r > r_|-, are : 



i?(r) = Cii?i(r) + C2i?2(r), for 



(2) 



i?i(r) = e°2" (r — r+) {r — r^)'^ HeunC(a, /3, 7, 5, 77, z) 
Rair) = e~ (r — r_|_)^'~ (r — 7'_)^HeunC(a,— /3, 7, S, ij, z), 



where z = - 

a = -2i I 



r— r-j- 



and the parameters are: 

^ ^ 2i{ui{a^+r.^}+am) 

-r_)uj,l3^ ^ 



1, 



7 = 



2i(ci; {a^ + r_-')+ani) 



1, 



[\-'iaujm-2uj^a^-2E 



-r 



4 (70; r_ — 270; r^+E 



^) r_ 7-^ — 4a^ (771+cj a) 
Here we have followed maple's internal rules when 
constructing the general solution of the differential equa- 
tion from the confluent Heun type. Accounting for the 
symmetries of the confluent Heun function, the solutions 

(3) coincide with those in [_ , ] (for uj replaced with ~uj). 

3 

The TRE has 3 singular points 7^. , r_|_, 00 and in order 
to fix the spectrum, one needs to impose specific bound- 
ary conditions on two of those singularities (i.e. to solve 
the central two-point connection problem [(if)]). Different 
boundary conditions on different pairs of singular points 
will specify different physics of the problem. In our case, 
we impose the black hole boundary conditions (BHBC) 
- waves going simultaneously into the event horizon (r+) 
and into infinity - following the same reasoning as in [Gf] 
where additional details can be found. Then, the BHBC 
read: 

1. BHBC on the KBH event horizon r+. 



For r — > r+, from r(t) = r+ + e ^("12) 
where 711,2 are the powers of the factors (r — r+)"i'2 
in R1.2, follows that for tti = 0, the only valid solu- 
tion in the whole interval (—00, 00) is R2, while for 
777, ^ 0, the solution R2 is valid for frequencies for 
which K(w) ^ {-JMTT^^)- This means that the ro- 
tation splits the area of validity of R2 into two and 
if this condition is not fulfilled, then the spectrum 
corresponds to waves going out of the horizon - a 
white hole. We won't pursue the spectrum in the 
case of a white hole, but it is important to keep in 



^ It is important to emphasize that the so obtained solutions can- 
not be used for extremal KBH (a = M) since in this case the 
differential equation is of the double confluent type and its treat- 
ment differs, so it is outside the scope of this work. 



mind that the boundary conditions correspond to 
a BH, only in the ranges of validity of each solu- 
tion. In our numerical work we use only R2 since 
the confluent Heun function in i?i is numerically 
unstable in maple. 

2. BHBC at infinity. 

At 7' — > 00, the solution is a linear combination 
of an ingoing (i?<-) and an outgoing (R^) wave: 
R = C<_ + C-i. R^ , where C<_ , C_j. are un- 
known constants and i?^ , R^ are found using the 
asymptotics of the confluent Heun function as de- 
fined in [27, 69]. 

To ensure only outgoing waves at infinity, one needs 
to have C<_ = 0. 

To achieve this, first one finds the direction of 
steepest descent in the complex plane for which 

lim = ^-'^^uM+2^-2^uJr ^ q ^^^^^ ^ero 

most quickly: sin(arg(a;) + arg(7')) = — f . This gives 
us a relation r =| r | e^/^*'^-* arg(<^) Q22]) between uj 
and r which is exact only if one uses the first term of 
the asymptotic series for the confluent Heun func- 
tion (i.e. HeunC ^ 1). More details about this 
approximation can be found in the next section. 

Then, it is enough to solve : 



2+iu>^£LSLl 

C^=r ' + 



HeunC(a,— /?, 7, 5, 77, z) = 0, (3) 



in order to completely specify the spectra 
{cj„,„i, with r = no e^/^"-* arg(a;) (^^^ 

\r\ — 110 as the actual numerical infinity and 
M =1/2). 



THE EPSILON-METHOD 

Equation (3) relies on the direction of steepest de- 
scent defined by the phase condition arg(r) + arg(tj) = 
3/27r. This approximate direction was chosen ignoring the 
higher terms in the asymptotic expansion of the solution 
around the infinity point, therefore, one can expect that 
the true path in the complex plane may not be a straight 
line but a curve. In principle, the spectrum should not 
depend on this curve, as long as r stays in the sector 
of the complex plane where lim = 0, i.e. as long as 

TT < arg(r)-|-arg(a;) < 27r, with only the convergence of the 
algorithm being affected. Numerical exploration of that 
limit evaluated with the first 3 terms in the asymptotic 
expansion of the appropriate confluent Heun functions for 
the modes w„, 77, = 0..18 when there is no BH rotation 
(a = 0), and for some modes when there is rotation, con- 
firms that indeed the limit remains approximately zero 
in the whole interval tt + 0.1 < arg(r)-|-arg(a;) < 27r — 0.1. 

The spectrum obtained numerically in this interval, 
however, depends in a nontrivial way on this curve. 
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The complications are partially due to the appearance 
of branch cuts in the numerical realization of the con- 
fluent Heun functions in MAPLE. The branching points 
of the confluent Heun function in the complex z-plane 
are found at the singular points z — 1 and z = oo. In 
MAPLE, the semi-infinite interval (1, oo) on the real axis 
is chosen as a branch cut. In the case of QNMs of non- 
rotating BHs [73, 74], it was observed that when those 
branch cuts are found near a frequency (since the radial 
variable, r depends on the frequency uj, the branch cuts 
appear also in the complex w-plane), they have serious 
effect on it, leading to the disappearance or translation 
of certain modes. As seen from the numerical study, this 
effect is likely due to a transition to another sheet of the 
multivalued function. 

As a way to find the correct sheet and remain on it, we 
introduced the epsilon-method which consists in adding 
a small variation (|e|<l) in the phase condition: 

arg(r)+arg(w)=:^-^^. (4) 

Using the e-method, one can change the location of the 
branch cut with respect to the eventual roots of the sys- 
tem and this way try to minimize the effect of the jump 
discontinuity of the radial function Using e, one can 
also explore the whole sector tt < arg(r)-|-arg(a;) < 2n 
, i.e. effectively moving r = |r|e"^''s('') jn the complex 
plane and this way test the numerical stability of the 
QNM spectrum, especially with respect to the position 
of the branch cuts of the radial function. 

Using the parameter e, the observed branch cuts in the 
realization of the confluent Heun function in MAPLE are 
as follows: 

1. For r~real. one encounters one of the branch cuts 
of the confluent Heun function. The equation 
of the line of this branch cut is: 3(w)/3fi(a;) = 
tan(3/27r -I- e7r/2) = — cot(e7r/2). This line rotates 
when e changes. 

2. If ~ 0, then one encounters the branch cut 
of the argument-function. In this case the branch 
cut is defined for = (—00,0). This branch 
cut, however, affects the solutions only very close 
to a = M where the frequencies can become almost 
real. 

3. If 3?(w) = and 5(a;) = 2n , n = 1,2,3.., then 
one can have 3(r) =0 for certain values of e and 
thus reach the branch cut of the confluent Heun 
function on the real axis. This condition can affect 
modes which are very near the imaginary axis (for 
example, similar condition holds around the alge- 
braically special mode for a nonrotating BH). 



* Here, the radial function refers to the solutions of the radial 
equation and not to the differential equation itself. 



4. Additional branch cuts may appear in the cases 
where 3(r'') = 0, for k-noninteger or complex 
(where z = 1 — ?■). Those branch cuts depend on the 
numerical realization of the confluent Heun func- 
tion in MAPLE and the equations of their lines can 
be obtained numerically from the values of the func- 
tion for each e. For example, one such branch cut 
was observed for 3(a;)/K(u;) « tan(1.419 + e7r/2). 

Although the study of the full effect from the move- 
ment of the branch cuts on the EM spectrum of KBH, 
as it was done in the case of gravitational perturbations 
of nonrotating black holes ([7.')]); is outside the scope of 
this work, some preliminary results on the issue can be 
found in the Appendix. 

NUMERICAL ALGORITHMS 

The spectral equations we need to solve to find the 
spectrum uJn,m{a) for M = 1/2 are Eqs.(l) and (3). This 
system represents a two-dimensional connected problem 
of two complex variables - the frequency uj and the sep- 
aration parameter E - and in both of its equations one 
encounters the confiuent Heun function and in the case 
of the TAE - their derivatives. 

A system like that cannot be easily solved by conven- 
tional methods like the Newton method and the Broy- 
den method, as outlined in [73] and [75], since they 
do not work well with the confluent Heun function in 
MAPLE. For this reason, our team developed a new 
method, namely the two-dimensional Miiller algorithm 
which proved to be much better adapted to work with 
those functions. The details of the algorithm can be 
found in [73-7-'')]; but for completeness, we will men- 
tion only that it relies on the Miiller method which is 
a quadratic generalization of the secant method having 
better convergence than the latter. The new algorithm 
does not need the evaluation of derivatives, thus avoid- 
ing one of the biggest problems when using the confluent 
Heun function in MAPLE. Clearly, in the system we solve 
the angular spectral equation (Eq. (1)) includes deriva- 
tives, but in this case, they remain in the domain \z\ < 1, 
where they can be evaluated correctly (for most values of 
the parameters) and with precision comparable to that 
of the radial function. It is important to note that both 
Ll),E are found directly from the spectral system (Eqs. 
[(1) and (3)]) and with equal precision. ^ ^ 



^ The algorithm is realized in maple code and the numbers pre- 
sented below are obtained using maple 13 on the computer clus- 
ter Physon. The software floating point number is set to 64 
(unless stated otherwise), the precision of the algorithm - to 15 
digits. 

® An important precaution when working with the confluent Heun 
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NUMERICAL RESULTS FOR 
ELECTROMAGNETIC QNMS 



While the evaluation of QNMs is not new to physics, 
the actual numbers published for EM perturbations 
of KBH are scarce. Because of this, for com- 
parison, we use the numbers published by Berti et 
al. [25, 56], the numerical data can be found on 
http://www.phy.olemiss.edu/~berti/qnms.html. Those 
numbers were obtained using the continued fraction 
method which is still considered as the most accurate 
method for obtaining the QNMs from the KBH. The 
available control frequencies are n = 0..6 for I ~ 1 and 
n = 0..3 for I = 2. Using those "control" numbers de- 
noted as w^^, one can easily check the precision of 
the method. 

The first 10 modes of the spectrum obtained using the 
new method in the interval a ~ [0,M) can be found 
on http://tcpa.uni-sofia.bg/conf/research. In the Ap- 
pendix, one can find some of the QNMs for specific values 
of a. 



Non-rotating BH 




It is already well known that when there is no rota- 
tion (a = 0) the electromagnetic QNMs come in pairs 
symmetrical to the imaginary axis aj„^„i — i|3^('^ri,m)| + 
iQ{(jJn,m) (n = 0,1.. numbering the mode). In this case 
the system reduces to one equation - the radial function 
(3) (for E = l{l + 1)J ^ 1,2..) solved here using the 
one-dimensional Miiller algorithm. 

The results can be seen on Fig. 1 a) and b), where 
we plotted the QNMs for m = 0,1,2 and I = 1,1 = 
2. The behavior of the modes resembles the behavior of 
the gravitational QNMs. A numerical comparison with 
the frequencies obtained by Berti et al. shows that the 
average deviation is — tOn.ml ~ 10~^° for the first 

4 modes (n = 0..3, I ~ 1,2). For modes with n > 3 
(i.e. n — 4. .6 for / = 1), however, there is an unexpected 
deviation which starts for n = 4 from \ujfj^ — uJ4^m\ ~ 
0.007 and grows to |cj,f,„ — wg,™! ~ 0.022 for the last 
available control mode. 

To study this systematic deviation, we employed the 
e-method to test the stability of those frequencies with 
respect to small deviations in the phase-condition. The 
results for e: 0, ±0.05, ±0.15 are plotted on Fig. 1 a) 



function in maple is that its precision or over-all behavior may 
depend on different factors which are not always under user's 
control [7;?]. From our observations, it seems that one can 
trust around 11-12 digits of the frequencies at the worst, usually 
around 13 digits. The points we present are the maximum we 
could get out of MAPLE, but future improvements in maple's code 
may significantly expand the area of application of the method 
and/or also its precision. 



Figure 1. QNMs for a = 0, for m = 0, 1, 2. The red diamonds 
are obtained for e = 0, the green crosses - for e — 0.05, the 
magenta diagonal crosses - for e = 0.15. With blue circles are 
the control frequencies uj^^m- Some points cannot be differed 
because for them, the numbers for different e coincide with 
precision higher than 10 digits 

and b). From there one can see that for n > 3 the best 
coincidence with the control frequencies a;^„j occurs for 
e = ±0.15, while for ti < 3 the modes obtained for the 
different values of e are equal and coincide with uj^ m- 

Similarly to the gravitational perturbation for nonro- 
tating BH ([73]), the dependence uj{e) in the electromag- 
netic case is not a trivial one. Here, only the case when 
there is no rotation (a = 0) allows studying in detail this 
dependence, since when there is rotation finding the roots 
of the system in a whole interval for e is strongly limited 
by the computational cost of the algorithm. Exploring 
a = 0, however, allow us to gather important intuition 
on the behavior of the QNMs under changes in e. The 
detailed study of this case can be found in the Appendix. 

The most important results from this study are: 

1. The modes u!m,n with n < N, where N depends 
on {to, e}, with high precision, do not depend on 
e and both signs of their real parts are obtainable 
with any sign of e. For modes with n > N, the rule 

is c^™,„(e)/™ = sgn(e)|5ft(w™,„)™| +i3(c^„,„)"\ 

2. The interval for e, where ujm,n with a certain sign 
of its real part can be found, shortens with the in- 
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crease of n. 

3. Studying ajm.„(e) for e G (—0.8. .0.8) shows that the 
effect due to e can be very small in certain ranges, 
but it is above the numerical error. 



Rotating KBH 

The results presented here were obtained for a = 
[0..M) for 3 different values of |e| = 0,0.05,0.15 (where 
the positive e are used for frequencies with positive real 
part and vice versa, negative e - for the frequencies with 
negative real parts). 

The results can be seen on the figures 2 — 11. 

When a ^ 0, the symmetry with respect to the imagi- 
nary axis aj^^„(0) = ±|5R(cJm.„)|-|-«3(ci;m,n) breaks down, 
but it is replaced by the symmetry: 



{K«^„), 3(£;,^„.„),m} ^ {-5ft«,„), -3(£;f„,„), -m}. 



where j — 1,2 coincides with the upper index of „(0). 
Thus, to study the complete behavior of the modes for 
a € [0, M), it is enough to trace both symmetric frequen- 
cies in the pair corresponding to each {m,n} for a = 0, 
for only m > (the index / here is omitted to simplify 
the notation, but everywhere in the text, if not explicitly 
stated otherwise, we compare only frequencies with the 
same I.) 

The parameter e, for rotating BH. has an even more 
significant role than the nonrotating case, since it does 
not merely translate modes with respect to each other, 
but for some modes, the frequencies obtained for the 3 
values of e have different behavior with respect to changes 
in the rotation. Some details can be found in the Ap- 
pendix. 

As in the case a ~ 0, one sees that for modes with n < 
^m,n(a)>w|^°«(a),w,°„^^(a)' coincide for equal {m,n} 
(Fig. 2), while when n > N, the modes „ obtained 
for different values of the parameter e differ (Fig. 3). TV 
depends on m, n and generally it is for N ~ 2. A when 
one of the modes (w^ „, wJJ;"^, w^^^) splits up from the 
rest (Note, here we discuss mostly the frequencies, but 
the separation parameters Em.n also depend on e as the 
figures show.). 

It is important to note that for m = 0, in the modes 
71 > 3, one observes loops. An example can be seen on 
Fig. 4. Those loops appear in all the higher modes, and 
their position depends on n. Because those loops require 
a finer structure of the plot (i.e. smaller step), on the 
plots Fig. 5, Fig. 6 and Fig. 9, we will plot only the 
points before the first loop observed in each curve. On 
Fig. 9 one can see all the results plotted together. 



0. 1 0.2 



(a)SR(L^0,3)(a) 



(b)3(L^0,3)(a) 
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(c)SR{£;o,3)(a) 



(d)3(£;o,3)(a) 



Figure 2. The plots depict the real and imaginary parts of 
cjo,3(a) and £"0,3 (a). Here, the red lines are the points ob- 
tained for 6 = 0, the blue crosses - those for e = 0.05 and the 
green diamonds - those for e = 0.15. The modes obtained for 
the 3 values of e coincide 




02 0.3 0.4 



(a)SR(LJo,4)(a) 



(b)3(Li;o,4)(a) 



J 




(c)SR(£;o,4)(a) 



(d)3(£;o,4)(a) 



Figure 3. The plots depict the real and imaginary parts 
ofa;o,4(a) and iJo,4(a). With red lines are denoted the points 
obtained for e = 0, with blue crosses - those for e = 0.05 and 
with green dashed lines - those for e = 0.15. The points for 
the different values of e differ 



where u)m,n{t) = n to avoid confusion with u)m,n{0') 



Figure 4. Example of the loops observed for m — 0. The 
figure shows the complex plots of tJo,3(i) and cjo,4(a), where 
the red lines are the points corresponding to e = 0, the green 
dashed line - those to e — 0.15. For n = 3, the results for 
e — 0,0.05,0.15 coincide and thus only the points for e = 
are plotted 
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Figure 5. Complex plots of a;o,,i(a) and Eo,n{a) for a = [0, M) 
the first 22 modes with both positive and negative real parts 



From the radial boundary conditions it follows that 
only frequencies for which "^{uS) ^ (0, —'m' 2Mr+ ) corre- 
spond to black hole boundary conditions. Figure 10 a) 
shows that the so obtained spectrum obeys this condi- 
tion. A deviation from this condition was observed in 
[61], where some of the frequencies describing primary 
jets crossed the line defined by ^"^ijUr^ ' thus corre- 
sponding to a white hole solution. For the QNM spec- 
trum, however, this is not the case and the spectrum 
corresponds to perturbation of a black hole. 

From the same figure one can see in the negative sector 
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Figure 6. Complex plots of cj().„(a) and Eo,n{a) for a — [0, M), 
1 = 2, e^O 




0,1 02 0,3 



Figure 7. The plots show the real and the imaginary parts of 
uji^„{a) and Ei,n{a) for a = [0, M) for the modes n = 6. .20. 
The black solid circle denotes a = 

of the plot that the real parts of the QNMs for increasing 
n seem to tend to the line —m^j^ — , which requires fur- 
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Figure 8. Complex plots of iOm.n{ci) and Em,n{a) for a = 
[0, M), n = 0..3. With blue is m = 0, with green m — 1 




Figure 9. A complex plot of all the uJm,n{a.) and Em,n{a) 
obtained for a = [0, M) for m = 0, 1, Z = 1, 2 n = 0..10 



ther investigation for n > 10. For the positive sector (i.e. 
the frequencies with positive real parts) we were not able 
to trace the frequencies with high n near a — M, thus 
we cannot confirm the relation 5R(a;) = m for a ^ M 
observed in [19]. 

Finally, obtaining the modes in the limit a k, M could 
be of serious interest, if one is to compare the EM QNMs 
with the spectra obtained from astrophysical objects, but 
it is also technically challenging. This happens because 
for a = M the TRE changes its type and near this limit 
the confluent Heun function becomes numerically unsta- 
ble since these functions are transforming to the double 
confluent Heun ones. Because of this, the examination of 




Figure 10. On the plots: with red lines - 5R(ciJi,„)(a) and 
9(cji,„)(a) for a = [0, Af), n = 0..10 and with blue dashed 
line — mTTTT — for m = 1 




Figure 11. On the plot ^{u)i^n){a) and 9(tJi_,i)(a) for a = 
[0.49, 0.4995] for the modes n = 0, 1, 2, m = 1 ' 

the limit a ^ M for modes with high n is impossible with 
current numerical realization of that function in MAPLE. 
For the lowest modes, however, the function is stable 
enough in the interval a £ [0.49, 0.4995] and the results 
of the numerical experiment for m = 1 are plotted on Fig. 
11. As expected, for n ~ 0, for a > 0.91 A/ the imaginary 
part of the frequency quickly tends to zero, thus proving 
that for extremal objects, the perturbations damp very 
slowly. The other two modes also seem to tend to zero, 
although somewhat slower than n = 0. In physical units, 
the difference between the 3 modes for a = 0.4995 is only 
6Hz (wi^i Rs 1.582kHz), but the damping times of the 
first mode is approximately 4.86 times bigger than that 
of the third and is i!^°™P ~ 4.2tos for KBH with mass 
A/ = 10A/q. The frequencies in physical units, for some 
other values of the rotational parameter, can be found in 
the tables I in the Appendix. 

While the analytical study of the extremal case is out- 
side the scope of this work, one can find such analytical 
treatment of the issue in [21, 24[. In both articles, one 
deals with approximations of the exact solutions of the 
radial equation, obtained under certain assumptions - 
be it through the continued fraction method in the limit 
a Af ([21]) or through a particular case in which the so- 
lution of the radial equation for a M can be written in 
terms of confluent hypergeometric functions ([24]). Both 
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methods seem to describe well the numerical results, but 
further investigation through the exact solutions written 
in terms of double confluent Heun functions and their 
properties, is needed to properly understand this inter- 
esting from observational point of view case. 

Algebraically special modes and branch cuts 

The algebraically special (AS) modes are obtained 
from the condition that the Starobinsky constant van- 
ishes ([1">]) and they correspond to the so called total 
transmission modes (TTM) ~ modes moving only in one 
direction: to the right or to the left. In the case of 
gravitational perturbations (s = —2) from nonrotating 
BH, since the 9*'' QNM coincides approximately with the 
theoretically expected purely imaginary AS mode, there 
were speculations that the two modes coincide (see [19] 
for a review, and also [76, 77]). A study of this mode in 
the case of gravitational perturbations of KBH showed 
numerical peculiarities as the "doublet" emerging from 
the "AS mode" for m > (see [19]) and also unexplained 
"spurious" modes, blamed to numerical inaccuracies. 

For electromagnetic perturbations, the possibility of 
appearance of the algebraically special modes has not 
been discussed much, since in the limit a — >^ 0, the 
Starobinsky constant does not vanish for purely imagi- 
nary modes (in fact, for a = 0, the Starobinsky constant 
does not depend on ut at all, see Eq. (60) [15] p. 392) and 
there appears to be no correlation between TTM and 
QNM modes [78]. There is, however, one important par- 
allel between the electromagnetic and the gravitational 
case. For the nonrotating gravitational case, Maassen 
van den Brink [76, 77] found that the peculiarities of the 
9*'' mode are due to the branch cut in the asymptotics of 
Regge- Wheeler potential, which the method of the con- 
tinued fraction is not adapted to handle. This result was 
confirmed by the use of the e-method in our previous 
work [74] where the AS character of the 9*'' mode was dis- 
proved. Using the e-method, one sees that this result is 
not limited to the gravitational case and the branch cuts 
play an important role for the electromagnetic QNMs as 
well. 

If one considers e as a parameter controlling the lo- 
cation of the branch cut with respect to certain QNM 
®, using the equations of the branch cuts discussed in 
section "The epsilon-method", one can find how "close" 
a certain branch cut is in the complex w-plane for each 
mode. In the gravitational case [74], the supposed AS 
mode is the one with the smallest real part and for it the 
value of e for which one observes the jump discontinuity 



Recall that r ~ exp(i arg(a;)) and thus by changing uJm.n and e, 
one changes the position in the r-complex plane 



is also very small. Therefore, one can expect that for 
this mode, very small variations in the phase-condition 
can change the leaf of the multivalued function and thus 
to lead to a different lo from the expected. 

In the EM case, one can also find a mode with a very 
small real part - n = 11 with 5R(a;o,ii) = .0215 (evaluated 
for e = 0.15), for which one encounters the jump discon- 
tinuity very close to the imaginary axis at e = 0.0024. 
In this case, one observes particularly interesting depen- 
dence on e - as showed on Fig. 1 a) - for e < 0.05 the 
mode n = 11 separates the lower QNM branch from the 
upper branch similarly to the way the so-called AS mode 
separates the QNM branches in the gravitational case 
[20], but if one uses e = 0.15 there is no such separation. 
Thus, for e = one finds a similarity between the EM 
and the gravitational cases. This similarity, however, is 
due to the appearance of branch cuts in the radial func- 
tion in both cases and not to some special properties of 
the mode in question (i.e. n = 11 for s = — 1 and n ~ 8 
for s ~ —2). This is because in the EM case, there is 
no theoretical justification for n = 11 to be an AS mode. 
In fact when the real part of that mode is very small 
(e < 0.05), its imaginary part deviates from the value of 
6i and vice versa - when the real part is not so small 
(e = 0.15), the imaginary part tends to 6i, see Fig. 1. 
Therefore, for all values of e, this mode deviates from a 
purely imaginary, integer number. 

An important note is that usually when studying 
purely imaginary frequencies, most methods do not work 
well, but in our case, there is no such problem, since we 
have a way to control some of the discontinuities in the 
solution of the radial function. If we are to continue the 
analogy with the gravitational case, studying how the 
mode 71 = 11 evolves with the increase of the rotation 
shows that it does no differ from the other modes. Al- 
though we didn't search specifically for doublets like the 
ones mentioned in [19], most of the modes with n > N 
can be considered as doublets, since we obtain distinct 
curves in the complex w— plane for different e. Further- 
more, in certain ranges of e we found two very close fre- 
quencies (i.e. corresponding to the same n) as roots of 
the system for the same e. The details for one of those 
cases can be found in the Appendix, but this also re- 
sembles the doublet found in the gravitational case, here 
found for n ~ 10. 

Focusing on the other peculiarity observed in [19] - 
the so called "spurious" modes the study of the fre- 
quencies in the interval e = —0.9. .0.9 for a ~ 
showed that indeed when varying e one may observe a 
number of frequencies around a certain mode. While 
the definition of "spurious" is unclear, in the case 
a = 0, n = 2, for example, one finds 3 frequen- 
cies (we omit here the separation constant which also 
differs): 0.3716378885 + 1.0257637188ii, 0.3610740790 + 
1.0392917852^,0.3495471352 + 1.0503751987i (corre- 
sponding to e = —0.22, —0.21, —0.2, respectively), where 
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the precision in all cases is 15 digits. Such examples oc- 
cur for all modes, with different roots found stable in 
different intervals for e. In some cases, one may consider 
those results as dependence of the frequency with e since 
the transition between those roots appears to be smooth 
^. In other cases, however, the transition appears to be 
step-wise in the scale of variation of e we studied. While 
there is no analytical explanation for this behavior yet, 
it appears to be related to the branch cut in the radial 
function which is moved by the parameter e. Thus, even 
if the phenomena of the "spurious" modes in [19] could 
indeed be blamed to numerical problems of the algorithm 
used in that article ([19]), exploring the roots with the 
e-method showed that a similar result can be obtained 
also in the EM case using a very different method. The 
reason for the observed peculiarities in the behavior of 
the spectrum of QNMs, due to the variation of epsilon, 
may be the complex character of the used analytical func- 
tions (the confluent Heun functions) in the vicinity of the 
irregular singular point r = oo in the complex r-plane. 

Considering all the numerical peculiarities demon- 
strated above, the use of the e-method poses a very 
serious question in front of the astrophysical applica- 
tion of that spectra - if one is to compare the numer- 
ical results with some observational frequencies, which 
e should be trusted? In our numerical experiments, we 
were able to obtain both the frequencies obtained with 
well-established methods with a precision bigger than 7 
digits, and also other, significantly deviating from them 
frequencies, which have qualitatively quite different be- 
havior with respect to changes in the rotation of the 
KBH. Both results are stable in different ranges for e, 
thus requiring new criteria for sifting out the physical 
modes based on better understanding of the behavior of 
the radial function in the complex plane of the radial 
variable. Such a study is outside the scope of the cur- 
rent work which aims to demonstrate the dependence of 
the method for obtaining the frequencies with respect to 
changes in the phase-condition and thus to provoke work 
in this area. 

CONCLUSION 

From the recent developments in the field of gravita- 
tional waves detection it is clear that finding the EM 
counterpart to those events can prove to be very useful. 
In this case, it is needed to better understand the fun- 
damental physics of quasi-normal ringing. In this paper. 



® In the mentioned above case (n = 2) one observes such smooth 
transitions: for e = -0.3. .-0.25 one finds the root .3901344442-1- 
.9754977543i which has 9 stable digits in this interval, which then 
goes smoothly to the "official" mode . 3495471352-1-1. 0503751987i, 
stable in the widest interval - e = — 0.2..0.7, through the numbers 
mentioned above for e = —0.22, —0.21, —0.2 



our team offered a new approach to finding the QNMs 
for the KBH based on directly solving the system ob- 
tained by the analytical solutions of the TRE and TAE 
in terms of the confluent Heun function. This approach 
has the advantage of being more traditional (i.e. impos- 
ing directly the corresponding boundary conditions on 
the exact analytical solutions of the problem) and hence 
it should allow better understanding of the peculiar prop- 
erties of the EM QNMs and the physics they imply. 

It was shown that using this approach, one can repro- 
duce the frequencies already obtained by other authors, 
but without relying on approximate methods. Partic- 
ularly important is the ability to impose the boundary 
condition directly on the solutions of the differential equa- 
tions. We require the standard regularity condition on 
the TAE and we explore in detail the radial boundary 
condition (the BHBC). Critical in it is the use of the di- 
rection of steepest descent, which secures the purely out- 
going wave at infinity. By using small deviations from 
this direction (and the phase-condition it defines), we 
were able to move around the branch cut in the solu- 
tions of the radial equation and thus to study its effect 
on the so obtained spectra. While this movement had 
no significant effect for the lower modes n < 3, for the 
higher modes it led to significant deviations from the al- 
ready published results. This behavior is persistent for 
the modes with m = 0, 1,2 and / = 1,2. This observa- 
tion raises the important question: What are the electro- 
magnetic QNMs for which one has to look in astrophys- 
ical data. Also interesting is that while the e-method 
leads to significant changes of the frequencies uim.m it 
affects much less the separation parameter Em,n which 
here for the first time was obtained directly as a solution 
of the two-dimensional system without any prior approx- 
imations for it. 

Another general result is that the confiuent Heun func- 
tion proved to be an effective tool for physical problems. 
Even though its maple realization still has many flaws, 
its precision proved to be good enough to repeat the al- 
ready published results, and also studying the solutions, 
we were able to reveal new properties of the numerical 
stability of the EM QNMs with respect to changes in the 
phase-condition. 

An interesting question is the results obtained using 
this method for a > M or the so called naked singular- 
ity regime. Preliminary results show that the method 
is applicable in this case as well and the results will be 
published elsewhere. 
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Tables of the obtained EM QNMs 

In the table I are presented some of the values obtained 
for the EM QNM, converted to physical units using the 
relations: 



phys 



2ttGM 



phys 



1 GM 



Note that in those formulas a factor of 2 is missing be- 
cause the EM QNMs were obtained for Mkbh = 1/2 
and not for Mkbh = 1- Then if M is the mass of 
the object in physical units, Mq - the mass of the Sun 
(M© = 1.98892 1030[fcg]) and G = 6.673 10-"[^], c = 
2.99792458 10^[to/s], one obtains 



phys 



32310 



phys 



0.4925 10-5 M/Mg 



The frequencies and the damping times in the table 
are calculated for M = IOMq. 



The e-method for a = 



Let us denote the dependence a;(e) as cj^^ (so that it 
differs from cj„(a)). 

First, one can compare the frequencies obtained for 
e = 0, 0.05, 0.15 for different m: 

for the case m = 0: \uj'i - a;° "5| « IQ-i^ fo^. n, 
but — a;°'^5| « 10^^*^ only for modes with n < 4, and 
1^0 05 -w0.i5| ^10-10 for n< 6. 

for the case m = 1: |w° - uj°J^\ « 10"^° and |a;° - 
a;0 i5| w 10-10 for n = 0..3 and n = 6. .11, and l^o o^ - 
a;° i5| w 10-12 for n < 12. 

for the case to 2: 1^° - « 10-^° for n = 0..3 

and n = 6.. 14, |w° - uj^-^^\ w 10-^° for n < 4. and 
1^0.05 _^o.i5|^ 10-10 forn< 6. 

The modes not enlisted above, such as and u>^^^ 
for n > A, m = 1, show significant deviation from the 
control results. 

Clearly, the modes for different m demonstrate differ- 
ent properties with respect to e. Such dependency on m 
is unexpected, since in equation (2), to is always coupled 
with a, so for a = 0, those frequencies should coincide. 
This indeed happens with precision lO-^^ for frequencies 
evaluated for the same e. For e 7^ 0, however, deviation 
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0.2 
0.6 


254.8509577191 
264.6648749590 
268.3651895175 


0.0275665662 
0.0277638854 
0.0294003836 


0.03 
0.31 
0.02 


254.8509577190 
334.8615860768 
610.5848397689 


0.0275665662 
0.0285332670 
0.0318459819 


0.34 
0.44 


254.8509577194 
251.7604458626 
291.2755338407 


0.0275665662 
0.0268474465 
0.0263079796 


0.03 
0.04 
0.06 


0.98 





















Table I. Table of the frequencies, cj'' in Hz , the damping times, t'' in milliseconds for some of the modes n — 0,4,7 and 
for some chosen values of the rotational parameter in the case I = 1. Presented is also S{ijj'^), the maximal difference between 
the modes obtained for the 3 values of e = 0, 0.05, 0.15 for each a, n, m. The numbers presented here correspond to 10 Mq. 
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Table II. Table of the separation parameter E for some of the modes n = 0, 4, 7 and for some chosen values of the rotational 
parameter in the case / = 1. Also S{E'^) is presented, the maximal difference between the E obtained for the 3 values of 
e = 0, 0.05, 0.15 for each a, n, m. 



may occur because different signs of e may lead to differ- 
ent roots. Moreover, one should not forget that the 



e = 0.15, where while most of the modes coincide with at least 
A peculiar case is the difference between m = and m = 1, for 11 digits, the frequencies n = 6. .11 show small deviation (r; 0.1), 
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modes are obtained through numerical algorithms with 
certain numerical zero. In this case, the factors in front 
of this number (i.e. m) may influence the final results 
for the frequencies, even though theoretically, one should 
expect the frequencies in the case a = not to depend on 
TO. Another hint for this is that the modes obtained from 
the radial spectral function for a = e = 0, 0.05, 0.15 co- 
incide with those obtained from the spectral system when 
a — > and that the points for a = seem not to deviate 
from the curves w„^„i(a). 

Another important result concerns which sign of 
|5R(a;„)| can be obtained using positive or negative e. 
The study of two cases (m = ±1) for n = 0..19 shows 
that while for lower modes n < N one can obtain both 
signs of |5R(w„)| for any e, for higher modes: w„(e)-''™ = 
sgn(e)|5R(w„)™| + i5(a;„)*" (where the index in and fin 
stay for initial and final value of the root-finding algo- 
rithm). It turns out that for n > iV, in order to find 
modes with positive real parts, one should use e > and 
vice versa. This result is also confirmed in the case a =/= 0. 
Here A'' depends on m and e and for to = ±1, N=3 or 
N=5. 

From the use of the e-method for gravitational pertur- 
bations (s = —2) of nonrotating BH, both in the RWE 
and in the TRE (a = 0) , it is known that deviations from 
the results obtained with the continued fraction method 
occur, because of branch cuts in the numerical realization 
of the confluent Heun function. Due to them one finds 
that ujn = ±s(7n(e)|3fi(a;„)| +iQ{u}n) is valid, in different 
intervals of e whose width depends on n. 

The results of a similar study of a;„(e) for the elec- 
tromagnetic QNMs, for a = 0, can be found on Fig. 12. 
The numerical data show that the change in uj„ due to 
the variation of e occurs in the 11*'' digit of the real and 
imaginary parts of uJn-, thus it is a very small effect. This 
effect, however, is two orders of magnitude above the nu- 
merical error for the confiuent Heun function expected to 
have at least 13 stable digits. The dependence cj„.„i(e) 
looks chaotic, but some kind of periodic behavior may be 
suspected. If one approximates a;„^,„(e) with a periodic 
function, its amplitude decreases with the increase of n 
(see Fig. 12 a), where the results for n = 0..5 are plotted). 
This seems to imply that the observed numerically ef- 
fect becomes less pronounced for the higher modes, even 
though one can expect an increase of the error of the nu- 
merical integration in the complex plane with n. This 
line of reasoning, however, applies only to the intervals 
where the dependence ojn (e) is smooth (see the discussion 
on page 12). 

Additional details on a;„^m(e) for n — 7 can be seen 
on Fig. 12 b). The observed behavior is not isolated. 



but it repeats in all modes and all to and also when one 
includes rotation (i.e. a > 0). 
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(a) 




which is probably due to numerical instability in the subroutine 
evaluating the confluent Heun function in maple. 



(b) 



Figure 12. a) Dashed line: the average deviation Aii(a;,i) = 
5 X 10^^|a;o „ — ^on^^l for n = 0..5, m = in the interval 
for e, where LOn^rn has at least 9 stable digits. Solid line: the 
width of that interval. Both the amplitude and the width of 
the intervals decrease with n. b) The dependence 5R((j;o,,i)(e) 
for m = 0,n = 7,a = for e = -0.28..0.18. The average 
deviation is 10~^^ 



The behavior of the modes for different e for a > 

Studying the modes n = 0..10 in the case of different 
{Z, TO, ?i} and different values of e one obtains the follow- 
ing results: 

• The case to = 0, Z = 1. For modes with n < 4, 
the results for the 3 values of e coincide and they 
can be seen on Fig. 13. Comparing with the 
control frequencies, one obtains ||a;„,o ~ '^j^oll 
lQ-^°,\\En,o - E^qW < 10-1" confirming that in 
this case the two methods - the continued fraction 
and our method- work comparably well. 

For n > 4, the frequencies obtained for differ- 



ent e split up in a way that uj (a) 



,0.0 
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0.25- % »„e 

0.25 0.30 0.35 0.40 0.45 0.50 2 2.05 2.10 

»(».) »(«.) 

Figure 13. Complex plots of a;o,n(a) and Eo,n{a) for the first 
4 modes n = 0..3, a — [0, M). The points obtained for e = 
0, 0.05, 0.15 coincide with more than 10 digits thus only one e 
is plotted. The black solid circle denotes a = 




O.lft O.IS 0.20 0.22 0.24 0.2(. 0.16 O.IS Q20 0.22 0.24 0.26 

»(».) »(«,) 

Figure 14. Complex plots of LOi,4{a) and £1,4(0) for a = 
[0, M), e = (red diamonds), e = 0.05(blue crosses) and 
e = 0.15 (green line). There is dramatical deviation of the 
points obtained for the different values of e 



while to (a) differs from them (excluded from this 
"rule" are 'n = 4, 5 for which w°(a) differs from the 
other two (Fig. 3). The numerical comparison with 
the control results w^„j, -B^^ show that they are 
closer to the results obtained for e = 0.15. 

Modes with 3fi(a;o,n) < 0: as mentioned in the previ- 
ous section, for modes with n > 2 the choice of sign 
of e becomes critical for the sign of the real part of 
the frequency and so we used e = 0, —0.05, —0.15. 

The results are symmetrical to those obtained 
for 3?(a;o,n) > with respect to the imaginary axis 
for uJo^n and with respect to the real axis for i?o,n 
(for 5R(wo,n) > 0, ^{Eo.n) < and vice versa). 

• The case to = 0, Z = 2. 

In this case, the modes obtained for e = 0, 0.05, 0.15 
do no split up until n = 7, where a;'' °^(a) = 
uj^-^^{a), while a;°(a) differs. This behavior contin- 
ues until n = 12, where it is a;°'^^(a) that deviates 
from the other two. The results for / = 2 can be 
seen on Fig. 6. 

• The case to = 1,1 = 1. 

The modes obtained for the three values of e coin- 
cide up to n < 4. For n = 4 (Fig. 14), the mode 
with e = differs. For n = 6, however, it is e = 0.15 
that differs, while the other two coincide and this 
behavior continues for higher modes. Note, here 
the deviation for different e is much more signifi- 
cant than the case to = (Fig. 14). 

These results show that the peculiarities observed 
when there is no rotation are inherited by the modes 



Using positive e in this case shows that for n = 3, the results for 
e = and e = 0.05 split up (e = 0.15 can not be used at ah), but 
it is for n = 6 where e = 0.05 can no longer be used for finding 
frequencies with negative real parts and only e = lead to the 
desired modes. 






/ 

0^ ■ : ■ , . , 

0.2 0.4 0.6 

Figure 15. Complex plot of the modes a;i,n(a) for m = 1,1 = 
l,n = 0..9 for two values of the parameter e: e = plotted 
with red lines and e = 0.15 plotted with green crosses. The 
red dashed line corresponds to the branch cut with equation 
9(a;)/K(a;) « tan(1.419 + e7r/2) for e = 0. One can see the 
effect of this branch cut on the modes with n > 3 



for a > 0. Studying the dependence of aj„j,„(e) for each 
mode for a certain interval of e is computationally ex- 
pensive, so we did it only for the case n = 10, a — 0.01. 
The results from the specific interval e ~ 0.0785. .0.088 
can be seen on Fig. 16. Besides the characteristic de- 
pendence uj{e), here, for the first time one finds two 
pairs of points [w„ , En] as roots for the same e (for 
e e [0.07862. .0.088034]), namely: 

[0.0680207667 -I- 5.1463791539j, 2.0021361645 -|- 0.0514150701*], 
[0.1419210235 -I- 5.0686957246i, 2.0028302584 -I- 0.0505980908i]. 

Such unexpected result has yet to be explained, consid- 
ering the big difference between the two frequencies. In 
any case it points to a behavior which must be studied 
more carefully in order to better understand the numeri- 
cal stability of the EM QNMs with respect to the branch 
cuts in the radial function. 
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0.080 0.082 0.084 0.086 0.088 



Figure 16. The dependence A(a;i,io)(e), where the solid hne 
denotes the real part, the dotted line - the imaginary part and 
A(tJi,io)(e) =a;|,io-^?;?o*^ The plot is for e = 0.0785..0.088, 
a = 0.01. The dependence of ii'i,io(e) is similar 
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